In [1]:

    

#Here we illustrate how to sample n-body particles with Galaxia.
"""
We will produce a simple n-body version of the LMC/SMC. We model it as a simple 3d Gaussian distribution.
Data that we put in:
SIMBAD:
smc_l = 302.8
smc_b = -44.3
lmc_l = 280.5
lmc_b = -32.9
GUMS paper:
smc_depth = 1.48 kpc
lmc_depth = 0.75 kpc
smc_diameter = 2.1kpc
lmc_diameter = 4.3 kpc
distance_lmc = 48.1 kpc
distance_smc = 60.6 kpc
mu_alpha_smc = 0.95 mas/yr
mu_delta_smc = -1.14 mas/yr
mu_alpha_lmc = 1.95 mas/yr
mu_delta_lmc = 0.43 mas/yr
vlos_smc = 158 km/s
vlos_lmc = 283 km/s
feh_smc = -1.2 +- 0.2 dex
feh_lmc = -0.75 +- 0.5
sfr_smc = constant
sfr_lmc = constant
mass_smc = (5.31 ± 0.05) × 10^8 Msun (stellar mass) https://ui.adsabs.harvard.edu/abs/2018MNRAS.478.5017R/abstract (total mass: 2.4 × 10^9 Msun)
mass_lmc = 5x10^9 Msun (same ratio as SMC total to stellar mass would be 5e9) total mass 1e11 Msun
# Scaling to GDR2 counts
mass_smc = (5.31 ± 0.05) × 10^8 Msun (stellar mass) https://ui.adsabs.harvard.edu/abs/2018MNRAS.478.5017R/abstract (total mass: 2.4 × 10^9 Msun)
mass_lmc = 2.5x10^9 Msun (we scale to GDR2 starcounts) (same ratio as SMC total to stellar mass would be 5e9) total mass 1e11 Msun

diameter_smc = 7000 lyr
diameter_lmc = 14000 lyr
velocity_dispersion_smc = 2 km/s guess

SMC = (x, y, z) in kpc
    (15.41505773, -36.45615675, -42.3529639)
 (v_x, v_y, v_z) in km / s
    (-30.27567249, -195.93860278, 132.71019931)>
LMC = (x, y, z) in kpc
    (-0.68927397, -39.70942142, -26.12549237)
 (v_x, v_y, v_z) in km / s
    (-79.46073554, -250.78236772, 205.32537167)>
"""


    

    
Out[1]:





'\nWe will produce a simple n-body version of the LMC/SMC. We model it as a simple 3d Gaussian distribution.\nData that we put in:\nSIMBAD:\nsmc_l = 302.8\nsmc_b = -44.3\nlmc_l = 280.5\nlmc_b = -32.9\nGUMS paper:\nsmc_depth = 1.48 kpc\nlmc_depth = 0.75 kpc\nsmc_diameter = 2.1kpc\nlmc_diameter = 4.3 kpc\ndistance_lmc = 48.1 kpc\ndistance_smc = 60.6 kpc\nmu_alpha_smc = 0.95 mas/yr\nmu_delta_smc = -1.14 mas/yr\nmu_alpha_lmc = 1.95 mas/yr\nmu_delta_lmc = 0.43 mas/yr\nvlos_smc = 158 km/s\nvlos_lmc = 283 km/s\nfeh_smc = -1.2 +- 0.2 dex\nfeh_lmc = -0.75 +- 0.5\nsfr_smc = constant\nsfr_lmc = constant\nmass_smc = (5.31 ± 0.05) × 10^8 Msun (stellar mass) https://ui.adsabs.harvard.edu/abs/2018MNRAS.478.5017R/abstract (total mass: 2.4 × 10^9 Msun)\nmass_lmc = 5x10^9 Msun (same ratio as SMC total to stellar mass would be 5e9) total mass 1e11 Msun\n# Scaling to GDR2 counts\nmass_smc = (5.31 ± 0.05) × 10^8 Msun (stellar mass) https://ui.adsabs.harvard.edu/abs/2018MNRAS.478.5017R/abstract (total mass: 2.4 × 10^9 Msun)\nmass_lmc = 2.5x10^9 Msun (we scale to GDR2 starcounts) (same ratio as SMC total to stellar mass would be 5e9) total mass 1e11 Msun\n\ndiameter_smc = 7000 lyr\ndiameter_lmc = 14000 lyr\nvelocity_dispersion_smc = 2 km/s guess\n\nSMC = (x, y, z) in kpc\n    (15.41505773, -36.45615675, -42.3529639)\n (v_x, v_y, v_z) in km / s\n    (-30.27567249, -195.93860278, 132.71019931)>\nLMC = (x, y, z) in kpc\n    (-0.68927397, -39.70942142, -26.12549237)\n (v_x, v_y, v_z) in km / s\n    (-79.46073554, -250.78236772, 205.32537167)>\n'

In [2]:

    

from astropy.coordinates import SkyCoord, ICRS, CartesianRepresentation, CartesianDifferential, Galactic, Galactocentric
import astropy.units as u


    

In [3]:

    

#SMC
c1 = SkyCoord(l = 302.8*u.degree, b = -44.3*u.degree, frame='galactic')
ra = c1.icrs.ra.value
dec = c1.icrs.dec.value
print("ra in deg:",ra)
print("dec in deg:",dec)
c = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, distance=60.6*u.kpc, frame='icrs',
             pm_ra_cosdec = 0.95*u.mas/u.yr, pm_dec = -1.14*u.mas/u.yr, radial_velocity = 158*u.km/u.s,
             galcen_distance = 8.0*u.kpc, z_sun = 15.0*u.pc,
             galcen_v_sun=CartesianDifferential(d_x=11.1*u.km/u.s, d_y=239.08*u.km/u.s, d_z=7.25*u.km/u.s))
c.galactocentric


    

    




ra in deg: 13.17927729718147
dec in deg: -72.82792474314047






    
Out[3]:





<SkyCoord (Galactocentric: galcen_coord=<ICRS Coordinate: (ra, dec) in deg
    (266.4051, -28.936175)>, galcen_distance=8.0 kpc, galcen_v_sun=(11.1, 239.08, 7.25) km / s, z_sun=15.0 pc, roll=0.0 deg): (x, y, z) in kpc
    (15.41505773, -36.45615675, -42.3529639)
 (v_x, v_y, v_z) in km / s
    (-30.27567249, -195.93860278, 132.71019931)>

In [4]:

    

#LMC
c1 = SkyCoord(l = 280.5*u.degree, b = -32.9*u.degree, frame='galactic')
ra = c1.icrs.ra.value
dec = c1.icrs.dec.value
print("ra in deg:",ra)
print("dec in deg:",dec)
c = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, distance=48.1*u.kpc, frame='icrs',
             pm_ra_cosdec = 1.95*u.mas/u.yr, pm_dec = 0.43*u.mas/u.yr, radial_velocity = 283*u.km/u.s,
             galcen_distance = 8.0*u.kpc, z_sun = 15.0*u.pc,
             galcen_v_sun=CartesianDifferential(d_x=11.1*u.km/u.s, d_y=239.08*u.km/u.s, d_z=7.25*u.km/u.s))
c.galactocentric


    

    




ra in deg: 80.8456130588062
dec in deg: -69.78267074987376






    
Out[4]:





<SkyCoord (Galactocentric: galcen_coord=<ICRS Coordinate: (ra, dec) in deg
    (266.4051, -28.936175)>, galcen_distance=8.0 kpc, galcen_v_sun=(11.1, 239.08, 7.25) km / s, z_sun=15.0 pc, roll=0.0 deg): (x, y, z) in kpc
    (-0.68927397, -39.70942142, -26.12549237)
 (v_x, v_y, v_z) in km / s
    (-79.46073554, -250.78236772, 205.32537167)>

In [5]:

    

%pylab inline
import ebf
import shutil
import subprocess
from astropy.io import fits
import os, sys
path = os.path.abspath('../library/')
if path not in sys.path:
    sys.path.append(path)
from convert_to_recarray import create_gdr2mock_mag_limited_survey_from_nbody


    

    




Populating the interactive namespace from numpy and matplotlib

In [6]:

    

# Need to specify where the GalaxiaData folder is and how to name the new simulation

nbody_folder = '/home/rybizki/Programme/GalaxiaData/'
folder = 'LSMC/'
nbody_filename = 'MCs'
folder_cat = '../output/MCs'


    

In [7]:

    

def create_n_body(nbody_folder,folder,filename):
    '''
    VERY IMPORTANT: NEEDS AN AGE-METALLICITY RELATION!
    '''
    # 2 input files for Galaxia need to be created
    folder_create = nbody_folder + 'nbody1/' + folder
    if os.path.exists(folder_create):
        shutil.rmtree(folder_create)
        os.mkdir(folder_create)
        print(folder_create, "existed and was recreated")
    else:
        os.mkdir(folder_create)

    # Age metallicity distribution is specified in this input file, using Model A from Just&Jahreiß 2010
    #fehdex = np.log10(zfe)
    # Here we can specify the AMR according to our needs
    stellar_masses = [6e8,3.8e9]
    particles_spawned = [10000,100000]
    feh_mean = [-1.2,-0.75]
    feh_disperson = [0.2,0.5]
    extends = [2.1,4.3]
    velocity_dispersion = [10,20]
    pos_x = [15.41505773,-0.68927397]
    pos_y = [-36.45615675,-39.70942142]
    pos_z = [-42.3529639,-26.12549237]
    vel_x = [-30.27567249,-79.46073554]
    vel_y = [-195.93860278,-250.78236772]
    vel_z = [132.71019931,205.32537167]
    names = ['smc','lmc']
    for j,name in enumerate(names):
        print(j,name)
        total_stellar_mass = stellar_masses[j]
        particles = particles_spawned[j]
        age = np.linspace(0,13.5,num = particles+1)[1:]
        fehdex = np.sort(np.random.normal(feh_mean[j],feh_disperson[j],particles))
        #ohdex = np.log10(zox)
        enhancement = np.zeros(len(age))
        mass = np.ones_like(age)*(total_stellar_mass/particles)
        age = age[::-1]
        depth = extends[j]
        vel_disp = velocity_dispersion[j]
        #################################

        # Here the ebf nbody input file is created
        pos = np.zeros((len(age),3))
        vel = np.zeros((len(age),3))
        for i in range(0,len(pos)):
            pos[i] = np.array([pos_x[j] + np.random.normal(0,depth/4), pos_y[j] + np.random.normal(0,depth/4), pos_z[j] + np.random.normal(0,depth/4)])
            vel[i] = np.array([vel_x[j] + np.random.normal(0,vel_disp), vel_y[j] + np.random.normal(0,vel_disp), vel_z[j] + np.random.normal(0,vel_disp)])
        print("median pos_x: ",np.median(pos[:,0]))
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '.ebf', '/mass', mass,'w')	
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '.ebf', '/feh', fehdex,'a')
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '.ebf', '/id', 1,'a')
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '.ebf', '/alpha', enhancement,'a')
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '.ebf', '/age', age,'a')
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '.ebf', '/pos3', pos,'a')
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '.ebf', '/vel3', vel,'a')

        # Preparing file for Enbid
        ps = np.concatenate((pos,vel),axis = 1)
        enbid_filename = '%s.dat' %(name)
        np.savetxt(enbid_filename,ps,fmt='%.6f')
        # Making the parameterfile for Enbid
        filedata = 'InitCondFile     %s\nICFormat                  0     \nSnapshotFileBase        _ph3\nSpatialScale            1 \nPartBoundary            7   \nNodeSplittingCriterion  1   \nCubicCells              1   \nMedianSplittingOn       0   \nTypeOfSmoothing      3\nDesNumNgb            64   \nVolCorr              1   \nTypeOfKernel           3 \nKernelBiasCorrection   1    \nAnisotropicKernel      0   \nAnisotropy             0   \nDesNumNgbA             128 \nTypeListOn        0\nPeriodicBoundaryOn 0 \n\n' %(enbid_filename)
        myparameterfile = "myparameterfile3"
        file = open(myparameterfile, "w")
        file.write(filedata)
        file.close()
        #Running Enbid
        args = ['./Enbid', myparameterfile]
        p = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE)	
        print("Enbid calculates smoothing length")
        (output, err) = p.communicate()
        # Writing to nbody smoothing length file 
        t = np.genfromtxt(enbid_filename + "_ph3.est",skip_header=1)
        d6 = np.zeros((len(pos),2))
        d6[:,0] = t[:,1]
        d6[:,1] = t[:,2]
        print(d6.shape,t.shape)
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '_d6n64_den.ebf', '/h_cubic', d6, 'w')

        # Same for 3d
        # Preparing file for Enbid
        enbid_filename = '%s3d.dat' %(name)
        np.savetxt(enbid_filename,pos,fmt='%.6f')
        # Making the parameterfile for Enbid
        filedata = 'InitCondFile     %s\nICFormat                  0     \nSnapshotFileBase        _ph3\nSpatialScale            1 \nPartBoundary            7   \nNodeSplittingCriterion  1   \nCubicCells              1   \nMedianSplittingOn       0   \nTypeOfSmoothing      3\nDesNumNgb            64   \nVolCorr              1   \nTypeOfKernel           3 \nKernelBiasCorrection   1    \nAnisotropicKernel      0   \nAnisotropy             0   \nDesNumNgbA             128 \nTypeListOn        0\nPeriodicBoundaryOn 0 \n\n' %(enbid_filename)
        myparameterfile = "myparameterfile3"
        file = open(myparameterfile, "w")
        file.write(filedata)
        file.close()
        #Running Enbid
        args = ['./Enbid3d', myparameterfile]
        p = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE)	
        print("Enbid calculates smoothing length")
        (output, err) = p.communicate()
        # Writing to nbody smoothing length file 
        t = np.genfromtxt(enbid_filename + "_ph3.est",skip_header=1)
        d3 = np.zeros((len(pos)))
        d3 = t[:,1]
        print(d3.shape,t.shape)
        ebf.write(nbody_folder + 'nbody1/' + folder + filename + name + '_d3n64_den.ebf', '/h_cubic', d3, 'w')

    # Here the file which tells Galaxia where to find the input file is created
    filedata = 'nbody1/%s\n         %d       1\n%s.ebf\n%s.ebf\n' %(folder,len(names), filename + names[0], filename + names[1])
    #filedata = 'nbody1/%s\n         %d       1\n%s.ebf\n' %(folder,1, filename + names[1])
    file = open(nbody_folder + "nbody1/filenames/" + filename + ".txt", "w")
    file.write(filedata)
    file.close()


    

In [ ]:

    

create_n_body(nbody_folder,folder,nbody_filename)


    

    




/home/rybizki/Programme/GalaxiaData/nbody1/LSMC/ existed and was recreated
0 smc
median pos_x:  15.417296816963976
Enbid calculates smoothing length
(10000, 2) (10000, 3)
Enbid calculates smoothing length
(10000,) (10000, 3)
1 lmc
median pos_x:  -0.6872082856228522
Enbid calculates smoothing length
(100000, 2) (100000, 3)
Enbid calculates smoothing length
(100000,) (100000, 3)

In [ ]:

    

for seed in np.arange(1):
    print(seed)
    create_gdr2mock_mag_limited_survey_from_nbody(nbody_filename = nbody_filename,nside = 512,
                                              outputDir = folder_cat + "_%d" %(seed),
                                   use_previous = False, delete_ebf = True,
                                  fSample = 0.1, make_likelihood_asessment=False, seed = seed)


    

    




0






    




/home/rybizki/anaconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88
  return f(*args, **kwds)
/home/rybizki/anaconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88
  return f(*args, **kwds)
/home/rybizki/anaconda3/lib/python3.6/site-packages/sklearn/ensemble/weight_boosting.py:29: DeprecationWarning: numpy.core.umath_tests is an internal NumPy module and should not be imported. It will be removed in a future NumPy release.
  from numpy.core.umath_tests import inner1d






    




../output/MCs_d/ existed and was recreated
Galaxia spawns catalogue

In [ ]:

    




    

In [ ]: